Single-cell RNA_Seq data analysis in R using Seurat

6 weeks HDM data vs 2 weeks HDM data and dog allergen data

JEBUNNAHAR MISHU

Loading the libraries

library(BiocManager)
library(devtools)
## Loading required package: usethis
## 
## Attaching package: 'devtools'
## The following object is masked from 'package:BiocManager':
## 
##     install
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 'SeuratObject' was built under R 4.4.0 but the current version is
## 4.4.1; it is recomended that you reinstall 'SeuratObject' as the ABI
## for R may have changed
## 
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
## 
##     intersect, t
library(SeuratObject)
library(remotes)
## 
## Attaching package: 'remotes'
## The following objects are masked from 'package:devtools':
## 
##     dev_package_deps, install_bioc, install_bitbucket, install_cran,
##     install_deps, install_dev, install_git, install_github,
##     install_gitlab, install_local, install_svn, install_url,
##     install_version, update_packages
## The following object is masked from 'package:usethis':
## 
##     git_credentials
library(targets)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(dplyr)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(Matrix)
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(pheatmap)
library(harmony)
## Loading required package: Rcpp

Read the file

base_path <- "/Users/cbm737/Documents/development"
# Set the directory
setwd(base_path)

Dirs <- list.dirs(path = base_path, full.names = FALSE, recursive = FALSE)
Dirs
## [1] "DA2WK_filtered_feature_bc_matrix"  "HDM2WK_filtered_feature_bc_matrix"
## [3] "HDM6WK_filtered_feature_bc_matrix"
#list.files(Dirs)

for (x in Dirs) {
  name <- gsub('_filtered_feature_bc_matrix','', x)
  
  # Construct file paths
  mtx_path <- file.path(base_path, x, "matrix.mtx.gz")
  features_path <- file.path(base_path, x, "features.tsv.gz")
  barcodes_path <- file.path(base_path, x, "barcodes.tsv.gz")
  
  # Read the data
  DAHDM_count <- ReadMtx(mtx = mtx_path, 
                         features = features_path, 
                         cells = barcodes_path)
  # create seurat object
  assign(name, CreateSeuratObject(count = DAHDM_count))
  
}

Merge datasets

HDMDA_Obj <- merge(HDM6WK, y = c(HDM2WK, DA2WK), add.cell.ids = c("HDM6_6WK", "HDM2_2WK", "DA2_2WK"), project = "HDMvsDA") 
#str(HDMDA_Obj)
#View(HDMDA_Obj@meta.data)
head(HDMDA_Obj@meta.data)
##                                orig.ident nCount_RNA nFeature_RNA
## HDM6_6WK_AAACCTGCAATTCCTT-1 SeuratProject       1859          962
## HDM6_6WK_AAACCTGCACATTTCT-1 SeuratProject       6871         2588
## HDM6_6WK_AAACCTGCATCCGGGT-1 SeuratProject       3226         1227
## HDM6_6WK_AAACCTGGTACAGTGG-1 SeuratProject       1637          842
## HDM6_6WK_AAACCTGGTATGCTTG-1 SeuratProject      20947         4600
## HDM6_6WK_AAACCTGGTCAGATAA-1 SeuratProject       1919         1185
HDMDA_Obj
## An object of class Seurat 
## 33502 features across 754038 samples within 1 assay 
## Active assay: RNA (33502 features, 0 variable features)
##  3 layers present: counts.1, counts.2, counts.3
# create the column
HDMDA_Obj$Sample <- rownames(HDMDA_Obj@meta.data)

# splite the column
HDMDA_Obj@meta.data <- separate(HDMDA_Obj@meta.data, col = 'Sample', into = c('Id', 'Time', 'Barcode'), sep = '_')
unique(HDMDA_Obj@meta.data$Id)
## [1] "HDM6" "HDM2" "DA2"
# Joinlayers
HDMDA_Obj <- JoinLayers(HDMDA_Obj)

Quality control

# percentage mitocondrial 
HDMDA_Obj[["PercentMT"]] <- PercentageFeatureSet(HDMDA_Obj, pattern = "^mt-")

# Removing Tcr  
HDMDA_Obj <- HDMDA_Obj[!grepl('^Tr[abdg][vjc]', rownames(HDMDA_Obj))] # mouse
HDMDA_Obj
## An object of class Seurat 
## 33248 features across 754038 samples within 1 assay 
## Active assay: RNA (33248 features, 0 variable features)
##  1 layer present: counts
head(HDMDA_Obj@meta.data)
##                                orig.ident nCount_RNA nFeature_RNA   Id Time
## HDM6_6WK_AAACCTGCAATTCCTT-1 SeuratProject       1859          962 HDM6  6WK
## HDM6_6WK_AAACCTGCACATTTCT-1 SeuratProject       6871         2588 HDM6  6WK
## HDM6_6WK_AAACCTGCATCCGGGT-1 SeuratProject       3226         1227 HDM6  6WK
## HDM6_6WK_AAACCTGGTACAGTGG-1 SeuratProject       1637          842 HDM6  6WK
## HDM6_6WK_AAACCTGGTATGCTTG-1 SeuratProject      20947         4600 HDM6  6WK
## HDM6_6WK_AAACCTGGTCAGATAA-1 SeuratProject       1919         1185 HDM6  6WK
##                                        Barcode PercentMT
## HDM6_6WK_AAACCTGCAATTCCTT-1 AAACCTGCAATTCCTT-1  2.850995
## HDM6_6WK_AAACCTGCACATTTCT-1 AAACCTGCACATTTCT-1  1.673701
## HDM6_6WK_AAACCTGCATCCGGGT-1 AAACCTGCATCCGGGT-1  4.463732
## HDM6_6WK_AAACCTGGTACAGTGG-1 AAACCTGGTACAGTGG-1  3.543067
## HDM6_6WK_AAACCTGGTATGCTTG-1 AAACCTGGTATGCTTG-1  1.446508
## HDM6_6WK_AAACCTGGTCAGATAA-1 AAACCTGGTCAGATAA-1  2.397082
# Visualize data before filtering as a violin plot
VlnPlot(HDMDA_Obj, features = c("nCount_RNA", "nFeature_RNA","PercentMT"), ncol = 3)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
## Warning: Removed 490493 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 490493 rows containing missing values or values outside the scale range
## (`geom_point()`).

# filtering
HDMDA_Obj_filtered <- subset(HDMDA_Obj, subset = nCount_RNA < 20000 & nFeature_RNA > 200 & PercentMT < 15)
HDMDA_Obj_filtered
## An object of class Seurat 
## 33248 features across 22709 samples within 1 assay 
## Active assay: RNA (33248 features, 0 variable features)
##  1 layer present: counts
# Visualize data before normalization as a violin plot
VlnPlot(HDMDA_Obj_filtered, features = c("nCount_RNA", "nFeature_RNA","PercentMT"), ncol = 3)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.

VlnPlot(HDMDA_Obj_filtered, features = c("nCount_RNA", "nFeature_RNA","PercentMT"), ncol = 3, pt.size = 0)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.

# Feature Scatter is used to visualize feature to feature relations
FeatureScatter(HDMDA_Obj_filtered, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") +
  geom_smooth(method = lm)
## `geom_smooth()` using formula = 'y ~ x'

Normalize the data, Find the variables features, Scale the data to avoid technical noise (batch effect) or biological sources ( different cell cycle), perform linear dimension reduction, Find neighbors, clusters and Run UMAP to visualization

HDMDA_Obj_Norm <- NormalizeData(object = HDMDA_Obj_filtered) %>%
                   FindVariableFeatures() %>% 
                    ScaleData() %>%
                     RunPCA() %>%
                      FindNeighbors(dims = 1:17) %>%
                       FindClusters(resolution = c(0.3, 0.5, 0.8)) %>%
                        RunUMAP(dims = 1:17)
## Normalizing layer: counts
## Finding variable features for layer counts
## Centering and scaling data matrix
## PC_ 1 
## Positive:  S100a10, Thy1, S100a11, AW112010, Crip1, Id2, S100a6, Lgals1, Cxcr6, Vim 
##     Maf, Tmsb10, Icos, S100a4, Rora, Tagln2, Emb, Ahnak, Rgs1, Il18r1 
##     Ramp1, Gapdh, Ifi27l2a, Hmgb2, Ly6a, Il7r, Ccr2, Lsp1, Ifngr1, Glrx 
## Negative:  H2-Aa, H2-Eb1, H2-Ab1, Cd74, Cd79a, H2-DMb2, Cd79b, Bank1, Ly6d, Fcmr 
##     Ms4a1, Ebf1, H2-Ob, Iglc3, Ly86, Mef2c, Iglc2, Cd83, Cd19, Blnk 
##     Apoe, Ifi30, Ciita, Ctsh, Fcrla, Ighd, H2-Oa, Igkc, Syk, Napsa 
## PC_ 2 
## Positive:  Tcf7, Klf2, Tmsb10, Rflnb, Ccr7, Cd79a, Igfbp4, Cd79b, S1pr1, Igkc 
##     Iglc2, Nsg2, Hmgn1, Ighm, Gimap6, H2-Ob, Slamf6, Iglc3, Ms4a1, Sell 
##     Lef1, Ebf1, Actn1, H2-DMb2, Ly6d, 1110034G24Rik, Bank1, Fcmr, Cd19, Mef2c 
## Negative:  Il1rl1, Fgl2, Lgals1, Ahnak, CAAA01147332.1, Dleu2, Atp5md, Tle5, Serinc3, Arl5a 
##     Rgs1, 2410006H16Rik, Aopep, Il13, Vim, Atp5mpl, Ndufb1-ps, Gas5, Lgmn, Neb 
##     Fabp5, Rbpj, Lgals3, Il2ra, Atp5o, Micos10, Cd200r1, Adam8, Castor1, Capg 
## PC_ 3 
## Positive:  Tle5, Il1rl1, Serinc3, Gm43305, Dleu2, Fgl2, 2410006H16Rik, 1810026B05Rik, Bank1, Cd79b 
##     Cd79a, Arl5a, Gas5, Fcmr, CAAA01147332.1, Plaat3, Neb, Tmem134, Ebf1, Ms4a1 
##     H2-DMb2, Castor1, Il13, Atp5md, Ly6d, Cd19, H2-Ob, Tut4, Atp5mpl, Il2ra 
## Negative:  Lyz2, Fcer1g, Tyrobp, Sirpa, Cd68, Mpeg1, Tgfbi, App, Cd63, Lrp1 
##     Fn1, Csf2rb, F10, Eps8, Cd9, Fcgr3, Cd300lf, Tnfaip2, Ctsk, Chil3 
##     Wfdc17, Cebpa, Alox5ap, Emilin2, Tgm2, Anpep, Gda, Rgl1, Clec4a3, Ccl6 
## PC_ 4 
## Positive:  Tle5, Gas5, Atp5md, Atp5mpl, 2410006H16Rik, Resf1, 1810026B05Rik, Ndufb1-ps, Atp5o, Il7r 
##     AY036118, Cct6a, Micos10, Rbis, Dleu2, Plaat3, Aopep, Snhg8, Micos13, Oip5os1 
##     Tmem134, Otulinl, Oga, Tut4, Rab5if, Tomm6, Shld1, Ciao2b, Gm2682, Trerf1 
## Negative:  Birc5, Pclaf, Ccna2, Ube2c, Nusap1, Cdca3, Spc24, Cdca8, Stmn1, Ccnb2 
##     Cks1b, Cenpf, Tpx2, Top2a, Esco2, Rrm2, Ska1, Kif11, Aurkb, Cenpe 
##     Cenpm, Cdk1, Mki67, Clspn, Hist1h2ab, Cit, Ckap2l, Mxd3, Knl1, Rad51ap1 
## PC_ 5 
## Positive:  Capg, S100a4, S100a11, S100a6, Vim, S100a10, Id2, Ramp1, Ucp2, Cd52 
##     Lgals1, Crip1, Cxcr6, Igfbp7, Rora, Aqp3, Glrx, Maf, Tagln2, Gpx1 
##     Fth1, Litaf, Ccr2, Lsp1, Ccr6, Pmm1, Socs2, Actg1, Rbpj, Tmem176b 
## Negative:  Tle5, Gas5, Atp5md, Atp5mpl, Tcf7, 2410006H16Rik, Resf1, Atp5o, 1810026B05Rik, Ndufb1-ps 
##     Satb1, Klf2, Gm2682, Cct6a, Klf3, Micos10, Plaat3, Ms4a4b, Rab5if, Otulinl 
##     Igfbp4, Rbis, AY036118, Oga, Snhg8, Micos13, Oip5os1, Ccr7, Shld1, Aopep
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 22709
## Number of edges: 758934
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9215
## Number of communities: 13
## Elapsed time: 2 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 22709
## Number of edges: 758934
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9010
## Number of communities: 19
## Elapsed time: 2 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 22709
## Number of edges: 758934
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8797
## Number of communities: 21
## Elapsed time: 2 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 16:30:06 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:30:06 Read 22709 rows and found 17 numeric columns
## 16:30:06 Using Annoy for neighbor search, n_neighbors = 30
## 16:30:06 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:30:08 Writing NN index file to temp file /var/folders/g4/2by5qwfx58b8zrlppt8n38300000gn/T//RtmpRVTgSd/file1bab75634c40
## 16:30:08 Searching Annoy index using 1 thread, search_k = 3000
## 16:30:11 Annoy recall = 100%
## 16:30:12 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 16:30:12 Initializing from normalized Laplacian + noise (using RSpectra)
## 16:30:13 Commencing optimization for 200 epochs, with 967326 positive edges
## 16:30:19 Optimization finished
str(HDMDA_Obj_Norm)
## Formal class 'Seurat' [package "SeuratObject"] with 13 slots
##   ..@ assays      :List of 1
##   .. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
##   .. .. .. ..@ layers    :List of 3
##   .. .. .. .. ..$ counts    :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. .. .. .. .. ..@ i       : int [1:38377842] 38 52 63 69 118 123 167 180 193 219 ...
##   .. .. .. .. .. .. ..@ p       : int [1:22710] 0 958 3543 4770 5611 6794 8671 10725 12629 14008 ...
##   .. .. .. .. .. .. ..@ Dim     : int [1:2] 33248 22709
##   .. .. .. .. .. .. ..@ Dimnames:List of 2
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. ..@ x       : num [1:38377842] 1 1 2 1 1 1 1 1 2 1 ...
##   .. .. .. .. .. .. ..@ factors : list()
##   .. .. .. .. ..$ data      :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. .. .. .. .. ..@ i       : int [1:38377842] 38 52 63 69 118 123 167 180 193 219 ...
##   .. .. .. .. .. .. ..@ p       : int [1:22710] 0 958 3543 4770 5611 6794 8671 10725 12629 14008 ...
##   .. .. .. .. .. .. ..@ Dim     : int [1:2] 33248 22709
##   .. .. .. .. .. .. ..@ Dimnames:List of 2
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. ..@ x       : num [1:38377842] 1.86 1.86 2.47 1.86 1.86 ...
##   .. .. .. .. .. .. ..@ factors : list()
##   .. .. .. .. ..$ scale.data: num [1:2000, 1:22709] -0.0124 -0.0119 -0.0201 -0.0132 -0.197 ...
##   .. .. .. ..@ cells     :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:22709, 1:3] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. .. .. ..$ dim     : int [1:2] 22709 3
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. ..@ features  :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:33248, 1:3] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:33248] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
##   .. .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. .. .. ..$ dim     : int [1:2] 33248 3
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:33248] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
##   .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. ..@ default   : int 1
##   .. .. .. ..@ assay.orig: chr(0) 
##   .. .. .. ..@ meta.data :'data.frame':  33248 obs. of  8 variables:
##   .. .. .. .. ..$ vf_vst_counts_mean                 : num [1:33248] 0.0 0.0 0.0 4.4e-05 0.0 ...
##   .. .. .. .. ..$ vf_vst_counts_variance             : num [1:33248] 0.0 0.0 0.0 4.4e-05 0.0 ...
##   .. .. .. .. ..$ vf_vst_counts_variance.expected    : num [1:33248] 0.0 0.0 0.0 4.4e-05 0.0 ...
##   .. .. .. .. ..$ vf_vst_counts_variance.standardized: num [1:33248] 0 0 0 1 0 ...
##   .. .. .. .. ..$ vf_vst_counts_variable             : logi [1:33248] FALSE FALSE FALSE FALSE FALSE TRUE ...
##   .. .. .. .. ..$ vf_vst_counts_rank                 : int [1:33248] NA NA NA NA NA 1194 NA NA NA NA ...
##   .. .. .. .. ..$ var.features                       : chr [1:33248] NA NA NA NA ...
##   .. .. .. .. ..$ var.features.rank                  : int [1:33248] NA NA NA NA NA 1194 NA NA NA NA ...
##   .. .. .. ..@ misc      : list()
##   .. .. .. ..@ key       : chr "rna_"
##   ..@ meta.data   :'data.frame': 22709 obs. of  11 variables:
##   .. ..$ orig.ident     : chr [1:22709] "SeuratProject" "SeuratProject" "SeuratProject" "SeuratProject" ...
##   .. ..$ nCount_RNA     : num [1:22709] 1859 6871 3226 1637 1919 ...
##   .. ..$ nFeature_RNA   : int [1:22709] 962 2588 1227 842 1185 1881 2058 1909 1385 1474 ...
##   .. ..$ Id             : chr [1:22709] "HDM6" "HDM6" "HDM6" "HDM6" ...
##   .. ..$ Time           : chr [1:22709] "6WK" "6WK" "6WK" "6WK" ...
##   .. ..$ Barcode        : chr [1:22709] "AAACCTGCAATTCCTT-1" "AAACCTGCACATTTCT-1" "AAACCTGCATCCGGGT-1" "AAACCTGGTACAGTGG-1" ...
##   .. ..$ PercentMT      : num [1:22709] 2.85 1.67 4.46 3.54 2.4 ...
##   .. ..$ RNA_snn_res.0.3: Factor w/ 13 levels "0","1","2","3",..: 5 5 9 1 4 4 4 5 5 9 ...
##   .. ..$ RNA_snn_res.0.5: Factor w/ 19 levels "0","1","2","3",..: 13 3 12 14 13 13 5 3 3 12 ...
##   .. ..$ RNA_snn_res.0.8: Factor w/ 21 levels "0","1","2","3",..: 16 3 14 15 16 13 5 3 3 14 ...
##   .. ..$ seurat_clusters: Factor w/ 21 levels "0","1","2","3",..: 16 3 14 15 16 13 5 3 3 14 ...
##   ..@ active.assay: chr "RNA"
##   ..@ active.ident: Factor w/ 21 levels "0","1","2","3",..: 16 3 14 15 16 13 5 3 3 14 ...
##   .. ..- attr(*, "names")= chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   ..@ graphs      :List of 2
##   .. ..$ RNA_nn :Formal class 'Graph' [package "SeuratObject"] with 7 slots
##   .. .. .. ..@ assay.used: chr "RNA"
##   .. .. .. ..@ i         : int [1:454180] 0 995 1024 1599 3985 4150 4274 4475 5376 17394 ...
##   .. .. .. ..@ p         : int [1:22710] 0 10 41 67 79 84 102 134 148 165 ...
##   .. .. .. ..@ Dim       : int [1:2] 22709 22709
##   .. .. .. ..@ Dimnames  :List of 2
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. ..@ x         : num [1:454180] 1 1 1 1 1 1 1 1 1 1 ...
##   .. .. .. ..@ factors   : list()
##   .. ..$ RNA_snn:Formal class 'Graph' [package "SeuratObject"] with 7 slots
##   .. .. .. ..@ assay.used: chr "RNA"
##   .. .. .. ..@ i         : int [1:1540577] 0 508 605 715 995 1024 1084 1156 1553 1575 ...
##   .. .. .. ..@ p         : int [1:22710] 0 42 111 220 251 290 366 465 522 602 ...
##   .. .. .. ..@ Dim       : int [1:2] 22709 22709
##   .. .. .. ..@ Dimnames  :List of 2
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. ..@ x         : num [1:1540577] 1 0.1111 0.0811 0.1111 0.25 ...
##   .. .. .. ..@ factors   : list()
##   ..@ neighbors   : list()
##   ..@ reductions  :List of 2
##   .. ..$ pca :Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
##   .. .. .. ..@ cell.embeddings           : num [1:22709, 1:50] 0.885 -1.466 -31.677 0.699 -1.37 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
##   .. .. .. ..@ feature.loadings          : num [1:2000, 1:50] -0.0236 -0.1313 -0.0414 0.0053 -0.1367 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:2000] "Chil3" "Cd74" "Lyz2" "Ccl5" ...
##   .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
##   .. .. .. ..@ feature.loadings.projected: num[0 , 0 ] 
##   .. .. .. ..@ assay.used                : chr "RNA"
##   .. .. .. ..@ global                    : logi FALSE
##   .. .. .. ..@ stdev                     : num [1:50] 6.34 5.34 4.72 4.46 4 ...
##   .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
##   .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ] 
##   .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ] 
##   .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ] 
##   .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ] 
##   .. .. .. ..@ misc                      :List of 1
##   .. .. .. .. ..$ total.variance: num 958
##   .. .. .. ..@ key                       : chr "PC_"
##   .. ..$ umap:Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
##   .. .. .. ..@ cell.embeddings           : num [1:22709, 1:2] -0.848 -4.571 -13.325 -0.866 -2.916 ...
##   .. .. .. .. ..- attr(*, "scaled:center")= num [1:2] 0.118 1.03
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. ..$ : chr [1:2] "umap_1" "umap_2"
##   .. .. .. ..@ feature.loadings          : num[0 , 0 ] 
##   .. .. .. ..@ feature.loadings.projected: num[0 , 0 ] 
##   .. .. .. ..@ assay.used                : chr "RNA"
##   .. .. .. ..@ global                    : logi TRUE
##   .. .. .. ..@ stdev                     : num(0) 
##   .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
##   .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ] 
##   .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ] 
##   .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ] 
##   .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ] 
##   .. .. .. ..@ misc                      : list()
##   .. .. .. ..@ key                       : chr "umap_"
##   ..@ images      : list()
##   ..@ project.name: chr "HDMvsDA"
##   ..@ misc        : list()
##   ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
##   .. ..$ : int [1:3] 5 0 2
##   ..@ commands    :List of 7
##   .. ..$ NormalizeData.RNA       :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "NormalizeData.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:29:42"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "NormalizeData(object = HDMDA_Obj_filtered)"
##   .. .. .. ..@ params     :List of 5
##   .. .. .. .. ..$ assay               : chr "RNA"
##   .. .. .. .. ..$ normalization.method: chr "LogNormalize"
##   .. .. .. .. ..$ scale.factor        : num 10000
##   .. .. .. .. ..$ margin              : num 1
##   .. .. .. .. ..$ verbose             : logi TRUE
##   .. ..$ FindVariableFeatures.RNA:Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "FindVariableFeatures.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:29:44"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "FindVariableFeatures(.)"
##   .. .. .. ..@ params     :List of 12
##   .. .. .. .. ..$ assay              : chr "RNA"
##   .. .. .. .. ..$ selection.method   : chr "vst"
##   .. .. .. .. ..$ loess.span         : num 0.3
##   .. .. .. .. ..$ clip.max           : chr "auto"
##   .. .. .. .. ..$ mean.function      :function (mat, display_progress)  
##   .. .. .. .. ..$ dispersion.function:function (mat, display_progress)  
##   .. .. .. .. ..$ num.bin            : num 20
##   .. .. .. .. ..$ binning.method     : chr "equal_width"
##   .. .. .. .. ..$ nfeatures          : num 2000
##   .. .. .. .. ..$ mean.cutoff        : num [1:2] 0.1 8
##   .. .. .. .. ..$ dispersion.cutoff  : num [1:2] 1 Inf
##   .. .. .. .. ..$ verbose            : logi TRUE
##   .. ..$ ScaleData.RNA           :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "ScaleData.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:29:44"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "ScaleData(.)"
##   .. .. .. ..@ params     :List of 9
##   .. .. .. .. ..$ assay             : chr "RNA"
##   .. .. .. .. ..$ model.use         : chr "linear"
##   .. .. .. .. ..$ use.umi           : logi FALSE
##   .. .. .. .. ..$ do.scale          : logi TRUE
##   .. .. .. .. ..$ do.center         : logi TRUE
##   .. .. .. .. ..$ scale.max         : num 10
##   .. .. .. .. ..$ block.size        : num 1000
##   .. .. .. .. ..$ min.cells.to.block: num 3000
##   .. .. .. .. ..$ verbose           : logi TRUE
##   .. ..$ RunPCA.RNA              :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "RunPCA.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:29:56"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "RunPCA(.)"
##   .. .. .. ..@ params     :List of 10
##   .. .. .. .. ..$ assay          : chr "RNA"
##   .. .. .. .. ..$ npcs           : num 50
##   .. .. .. .. ..$ rev.pca        : logi FALSE
##   .. .. .. .. ..$ weight.by.var  : logi TRUE
##   .. .. .. .. ..$ verbose        : logi TRUE
##   .. .. .. .. ..$ ndims.print    : int [1:5] 1 2 3 4 5
##   .. .. .. .. ..$ nfeatures.print: num 30
##   .. .. .. .. ..$ reduction.name : chr "pca"
##   .. .. .. .. ..$ reduction.key  : chr "PC_"
##   .. .. .. .. ..$ seed.use       : num 42
##   .. ..$ FindNeighbors.RNA.pca   :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "FindNeighbors.RNA.pca"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:29:58"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "FindNeighbors(., dims = 1:17)"
##   .. .. .. ..@ params     :List of 16
##   .. .. .. .. ..$ reduction      : chr "pca"
##   .. .. .. .. ..$ dims           : int [1:17] 1 2 3 4 5 6 7 8 9 10 ...
##   .. .. .. .. ..$ assay          : chr "RNA"
##   .. .. .. .. ..$ k.param        : num 20
##   .. .. .. .. ..$ return.neighbor: logi FALSE
##   .. .. .. .. ..$ compute.SNN    : logi TRUE
##   .. .. .. .. ..$ prune.SNN      : num 0.0667
##   .. .. .. .. ..$ nn.method      : chr "annoy"
##   .. .. .. .. ..$ n.trees        : num 50
##   .. .. .. .. ..$ annoy.metric   : chr "euclidean"
##   .. .. .. .. ..$ nn.eps         : num 0
##   .. .. .. .. ..$ verbose        : logi TRUE
##   .. .. .. .. ..$ do.plot        : logi FALSE
##   .. .. .. .. ..$ graph.name     : chr [1:2] "RNA_nn" "RNA_snn"
##   .. .. .. .. ..$ l2.norm        : logi FALSE
##   .. .. .. .. ..$ cache.index    : logi FALSE
##   .. ..$ FindClusters            :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "FindClusters"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:06"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "FindClusters(., resolution = c(0.3, 0.5, 0.8))"
##   .. .. .. ..@ params     :List of 11
##   .. .. .. .. ..$ graph.name      : chr "RNA_snn"
##   .. .. .. .. ..$ cluster.name    : chr [1:3] "RNA_snn_res.0.3" "RNA_snn_res.0.5" "RNA_snn_res.0.8"
##   .. .. .. .. ..$ modularity.fxn  : num 1
##   .. .. .. .. ..$ resolution      : num [1:3] 0.3 0.5 0.8
##   .. .. .. .. ..$ method          : chr "matrix"
##   .. .. .. .. ..$ algorithm       : num 1
##   .. .. .. .. ..$ n.start         : num 10
##   .. .. .. .. ..$ n.iter          : num 10
##   .. .. .. .. ..$ random.seed     : num 0
##   .. .. .. .. ..$ group.singletons: logi TRUE
##   .. .. .. .. ..$ verbose         : logi TRUE
##   .. ..$ RunUMAP.RNA.pca         :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "RunUMAP.RNA.pca"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:19"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "RunUMAP(., dims = 1:17)"
##   .. .. .. ..@ params     :List of 25
##   .. .. .. .. ..$ dims                : int [1:17] 1 2 3 4 5 6 7 8 9 10 ...
##   .. .. .. .. ..$ reduction           : chr "pca"
##   .. .. .. .. ..$ assay               : chr "RNA"
##   .. .. .. .. ..$ slot                : chr "data"
##   .. .. .. .. ..$ umap.method         : chr "uwot"
##   .. .. .. .. ..$ return.model        : logi FALSE
##   .. .. .. .. ..$ n.neighbors         : int 30
##   .. .. .. .. ..$ n.components        : int 2
##   .. .. .. .. ..$ metric              : chr "cosine"
##   .. .. .. .. ..$ learning.rate       : num 1
##   .. .. .. .. ..$ min.dist            : num 0.3
##   .. .. .. .. ..$ spread              : num 1
##   .. .. .. .. ..$ set.op.mix.ratio    : num 1
##   .. .. .. .. ..$ local.connectivity  : int 1
##   .. .. .. .. ..$ repulsion.strength  : num 1
##   .. .. .. .. ..$ negative.sample.rate: int 5
##   .. .. .. .. ..$ uwot.sgd            : logi FALSE
##   .. .. .. .. ..$ seed.use            : int 42
##   .. .. .. .. ..$ angular.rp.forest   : logi FALSE
##   .. .. .. .. ..$ densmap             : logi FALSE
##   .. .. .. .. ..$ dens.lambda         : num 2
##   .. .. .. .. ..$ dens.frac           : num 0.3
##   .. .. .. .. ..$ dens.var.shift      : num 0.1
##   .. .. .. .. ..$ verbose             : logi TRUE
##   .. .. .. .. ..$ reduction.name      : chr "umap"
##   ..@ tools       : list()
head(HDMDA_Obj_Norm@meta.data)
##                                orig.ident nCount_RNA nFeature_RNA   Id Time
## HDM6_6WK_AAACCTGCAATTCCTT-1 SeuratProject       1859          962 HDM6  6WK
## HDM6_6WK_AAACCTGCACATTTCT-1 SeuratProject       6871         2588 HDM6  6WK
## HDM6_6WK_AAACCTGCATCCGGGT-1 SeuratProject       3226         1227 HDM6  6WK
## HDM6_6WK_AAACCTGGTACAGTGG-1 SeuratProject       1637          842 HDM6  6WK
## HDM6_6WK_AAACCTGGTCAGATAA-1 SeuratProject       1919         1185 HDM6  6WK
## HDM6_6WK_AAACCTGGTCCAGTAT-1 SeuratProject       4150         1881 HDM6  6WK
##                                        Barcode PercentMT RNA_snn_res.0.3
## HDM6_6WK_AAACCTGCAATTCCTT-1 AAACCTGCAATTCCTT-1  2.850995               4
## HDM6_6WK_AAACCTGCACATTTCT-1 AAACCTGCACATTTCT-1  1.673701               4
## HDM6_6WK_AAACCTGCATCCGGGT-1 AAACCTGCATCCGGGT-1  4.463732               8
## HDM6_6WK_AAACCTGGTACAGTGG-1 AAACCTGGTACAGTGG-1  3.543067               0
## HDM6_6WK_AAACCTGGTCAGATAA-1 AAACCTGGTCAGATAA-1  2.397082               3
## HDM6_6WK_AAACCTGGTCCAGTAT-1 AAACCTGGTCCAGTAT-1  2.024096               3
##                             RNA_snn_res.0.5 RNA_snn_res.0.8 seurat_clusters
## HDM6_6WK_AAACCTGCAATTCCTT-1              12              15              15
## HDM6_6WK_AAACCTGCACATTTCT-1               2               2               2
## HDM6_6WK_AAACCTGCATCCGGGT-1              11              13              13
## HDM6_6WK_AAACCTGGTACAGTGG-1              13              14              14
## HDM6_6WK_AAACCTGGTCAGATAA-1              12              15              15
## HDM6_6WK_AAACCTGGTCCAGTAT-1              12              12              12
ElbowPlot(HDMDA_Obj_Norm)

PCHeatmap(HDMDA_Obj_Norm, dim = 1:6, cells = 500, balanced =TRUE, ncol = 3)

DimPlot(HDMDA_Obj_Norm, reduction = "umap", group.by = "RNA_snn_res.0.5", label = TRUE, label.box = TRUE)+ labs(title = "Without Integration")

DimPlot(HDMDA_Obj_Norm, reduction = "umap", group.by = "RNA_snn_res.0.5", split.by = "Id", label = TRUE, label.box = TRUE)+ labs(title = "Without Integration")

Normalize the data, Find the variables features, Scale the data to avoid technical noise (batch effect) or biological sources ( different cell cycle), Perform linear dimension reduction, Integration using Harmony to avoid technical noise (batch effect), Find neighbors, clusters and Run UMAP to visualization

HDMDA_Obj_Harm <- NormalizeData(object = HDMDA_Obj_filtered) %>%
                   FindVariableFeatures() %>% 
                    ScaleData() %>%
                     RunPCA() %>%
                      RunHarmony(group.by.vars = "Id") %>%
                       FindNeighbors(reduction = "harmony", dims = 1:17) %>%
                        FindClusters(resolution = 0.5) %>%
                         RunUMAP(reduction = "harmony", dims = 1:17)
## Normalizing layer: counts
## Finding variable features for layer counts
## Centering and scaling data matrix
## PC_ 1 
## Positive:  S100a10, Thy1, S100a11, AW112010, Crip1, Id2, S100a6, Lgals1, Cxcr6, Vim 
##     Maf, Tmsb10, Icos, S100a4, Rora, Tagln2, Emb, Ahnak, Rgs1, Il18r1 
##     Ramp1, Gapdh, Ifi27l2a, Hmgb2, Ly6a, Il7r, Ccr2, Lsp1, Ifngr1, Glrx 
## Negative:  H2-Aa, H2-Eb1, H2-Ab1, Cd74, Cd79a, H2-DMb2, Cd79b, Bank1, Ly6d, Fcmr 
##     Ms4a1, Ebf1, H2-Ob, Iglc3, Ly86, Mef2c, Iglc2, Cd83, Cd19, Blnk 
##     Apoe, Ifi30, Ciita, Ctsh, Fcrla, Ighd, H2-Oa, Igkc, Syk, Napsa 
## PC_ 2 
## Positive:  Tcf7, Klf2, Tmsb10, Rflnb, Ccr7, Cd79a, Igfbp4, Cd79b, S1pr1, Igkc 
##     Iglc2, Nsg2, Hmgn1, Ighm, Gimap6, H2-Ob, Slamf6, Iglc3, Ms4a1, Sell 
##     Lef1, Ebf1, Actn1, H2-DMb2, Ly6d, 1110034G24Rik, Bank1, Fcmr, Cd19, Mef2c 
## Negative:  Il1rl1, Fgl2, Lgals1, Ahnak, CAAA01147332.1, Dleu2, Atp5md, Tle5, Serinc3, Arl5a 
##     Rgs1, 2410006H16Rik, Aopep, Il13, Vim, Atp5mpl, Ndufb1-ps, Gas5, Lgmn, Neb 
##     Fabp5, Rbpj, Lgals3, Il2ra, Atp5o, Micos10, Cd200r1, Adam8, Castor1, Capg 
## PC_ 3 
## Positive:  Tle5, Il1rl1, Serinc3, Gm43305, Dleu2, Fgl2, 2410006H16Rik, 1810026B05Rik, Bank1, Cd79b 
##     Cd79a, Arl5a, Gas5, Fcmr, CAAA01147332.1, Plaat3, Neb, Tmem134, Ebf1, Ms4a1 
##     H2-DMb2, Castor1, Il13, Atp5md, Ly6d, Cd19, H2-Ob, Tut4, Atp5mpl, Il2ra 
## Negative:  Lyz2, Fcer1g, Tyrobp, Sirpa, Cd68, Mpeg1, Tgfbi, App, Cd63, Lrp1 
##     Fn1, Csf2rb, F10, Eps8, Cd9, Fcgr3, Cd300lf, Tnfaip2, Ctsk, Chil3 
##     Wfdc17, Cebpa, Alox5ap, Emilin2, Tgm2, Anpep, Gda, Rgl1, Clec4a3, Ccl6 
## PC_ 4 
## Positive:  Tle5, Gas5, Atp5md, Atp5mpl, 2410006H16Rik, Resf1, 1810026B05Rik, Ndufb1-ps, Atp5o, Il7r 
##     AY036118, Cct6a, Micos10, Rbis, Dleu2, Plaat3, Aopep, Snhg8, Micos13, Oip5os1 
##     Tmem134, Otulinl, Oga, Tut4, Rab5if, Tomm6, Shld1, Ciao2b, Gm2682, Trerf1 
## Negative:  Birc5, Pclaf, Ccna2, Ube2c, Nusap1, Cdca3, Spc24, Cdca8, Stmn1, Ccnb2 
##     Cks1b, Cenpf, Tpx2, Top2a, Esco2, Rrm2, Ska1, Kif11, Aurkb, Cenpe 
##     Cenpm, Cdk1, Mki67, Clspn, Hist1h2ab, Cit, Ckap2l, Mxd3, Knl1, Rad51ap1 
## PC_ 5 
## Positive:  Capg, S100a4, S100a11, S100a6, Vim, S100a10, Id2, Ramp1, Ucp2, Cd52 
##     Lgals1, Crip1, Cxcr6, Igfbp7, Rora, Aqp3, Glrx, Maf, Tagln2, Gpx1 
##     Fth1, Litaf, Ccr2, Lsp1, Ccr6, Pmm1, Socs2, Actg1, Rbpj, Tmem176b 
## Negative:  Tle5, Gas5, Atp5md, Atp5mpl, Tcf7, 2410006H16Rik, Resf1, Atp5o, 1810026B05Rik, Ndufb1-ps 
##     Satb1, Klf2, Gm2682, Cct6a, Klf3, Micos10, Plaat3, Ms4a4b, Rab5if, Otulinl 
##     Igfbp4, Rbis, AY036118, Oga, Snhg8, Micos13, Oip5os1, Ccr7, Shld1, Aopep
## Transposing data matrix
## Initializing state using k-means centroids initialization
## Harmony 1/10
## Harmony 2/10
## Harmony converged after 2 iterations
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 22709
## Number of edges: 743589
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8794
## Number of communities: 16
## Elapsed time: 2 seconds
## 16:30:51 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:30:51 Read 22709 rows and found 17 numeric columns
## 16:30:51 Using Annoy for neighbor search, n_neighbors = 30
## 16:30:51 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:30:51 Writing NN index file to temp file /var/folders/g4/2by5qwfx58b8zrlppt8n38300000gn/T//RtmpRVTgSd/file1bab12d98dec
## 16:30:51 Searching Annoy index using 1 thread, search_k = 3000
## 16:30:55 Annoy recall = 100%
## 16:30:56 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 16:30:56 Initializing from normalized Laplacian + noise (using RSpectra)
## 16:30:57 Commencing optimization for 200 epochs, with 983656 positive edges
## 16:31:03 Optimization finished
str(HDMDA_Obj_Harm)
## Formal class 'Seurat' [package "SeuratObject"] with 13 slots
##   ..@ assays      :List of 1
##   .. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
##   .. .. .. ..@ layers    :List of 3
##   .. .. .. .. ..$ counts    :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. .. .. .. .. ..@ i       : int [1:38377842] 38 52 63 69 118 123 167 180 193 219 ...
##   .. .. .. .. .. .. ..@ p       : int [1:22710] 0 958 3543 4770 5611 6794 8671 10725 12629 14008 ...
##   .. .. .. .. .. .. ..@ Dim     : int [1:2] 33248 22709
##   .. .. .. .. .. .. ..@ Dimnames:List of 2
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. ..@ x       : num [1:38377842] 1 1 2 1 1 1 1 1 2 1 ...
##   .. .. .. .. .. .. ..@ factors : list()
##   .. .. .. .. ..$ data      :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. .. .. .. .. ..@ i       : int [1:38377842] 38 52 63 69 118 123 167 180 193 219 ...
##   .. .. .. .. .. .. ..@ p       : int [1:22710] 0 958 3543 4770 5611 6794 8671 10725 12629 14008 ...
##   .. .. .. .. .. .. ..@ Dim     : int [1:2] 33248 22709
##   .. .. .. .. .. .. ..@ Dimnames:List of 2
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. ..@ x       : num [1:38377842] 1.86 1.86 2.47 1.86 1.86 ...
##   .. .. .. .. .. .. ..@ factors : list()
##   .. .. .. .. ..$ scale.data: num [1:2000, 1:22709] -0.0124 -0.0119 -0.0201 -0.0132 -0.197 ...
##   .. .. .. ..@ cells     :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:22709, 1:3] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. .. .. ..$ dim     : int [1:2] 22709 3
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. ..@ features  :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:33248, 1:3] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:33248] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
##   .. .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. .. .. ..$ dim     : int [1:2] 33248 3
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:33248] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
##   .. .. .. .. .. .. ..$ : chr [1:3] "counts" "data" "scale.data"
##   .. .. .. ..@ default   : int 1
##   .. .. .. ..@ assay.orig: chr(0) 
##   .. .. .. ..@ meta.data :'data.frame':  33248 obs. of  8 variables:
##   .. .. .. .. ..$ vf_vst_counts_mean                 : num [1:33248] 0.0 0.0 0.0 4.4e-05 0.0 ...
##   .. .. .. .. ..$ vf_vst_counts_variance             : num [1:33248] 0.0 0.0 0.0 4.4e-05 0.0 ...
##   .. .. .. .. ..$ vf_vst_counts_variance.expected    : num [1:33248] 0.0 0.0 0.0 4.4e-05 0.0 ...
##   .. .. .. .. ..$ vf_vst_counts_variance.standardized: num [1:33248] 0 0 0 1 0 ...
##   .. .. .. .. ..$ vf_vst_counts_variable             : logi [1:33248] FALSE FALSE FALSE FALSE FALSE TRUE ...
##   .. .. .. .. ..$ vf_vst_counts_rank                 : int [1:33248] NA NA NA NA NA 1194 NA NA NA NA ...
##   .. .. .. .. ..$ var.features                       : chr [1:33248] NA NA NA NA ...
##   .. .. .. .. ..$ var.features.rank                  : int [1:33248] NA NA NA NA NA 1194 NA NA NA NA ...
##   .. .. .. ..@ misc      : list()
##   .. .. .. ..@ key       : chr "rna_"
##   ..@ meta.data   :'data.frame': 22709 obs. of  9 variables:
##   .. ..$ orig.ident     : chr [1:22709] "SeuratProject" "SeuratProject" "SeuratProject" "SeuratProject" ...
##   .. ..$ nCount_RNA     : num [1:22709] 1859 6871 3226 1637 1919 ...
##   .. ..$ nFeature_RNA   : int [1:22709] 962 2588 1227 842 1185 1881 2058 1909 1385 1474 ...
##   .. ..$ Id             : chr [1:22709] "HDM6" "HDM6" "HDM6" "HDM6" ...
##   .. ..$ Time           : chr [1:22709] "6WK" "6WK" "6WK" "6WK" ...
##   .. ..$ Barcode        : chr [1:22709] "AAACCTGCAATTCCTT-1" "AAACCTGCACATTTCT-1" "AAACCTGCATCCGGGT-1" "AAACCTGGTACAGTGG-1" ...
##   .. ..$ PercentMT      : num [1:22709] 2.85 1.67 4.46 3.54 2.4 ...
##   .. ..$ RNA_snn_res.0.5: Factor w/ 16 levels "0","1","2","3",..: 6 1 11 5 6 7 2 4 1 11 ...
##   .. ..$ seurat_clusters: Factor w/ 16 levels "0","1","2","3",..: 6 1 11 5 6 7 2 4 1 11 ...
##   ..@ active.assay: chr "RNA"
##   ..@ active.ident: Factor w/ 16 levels "0","1","2","3",..: 6 1 11 5 6 7 2 4 1 11 ...
##   .. ..- attr(*, "names")= chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   ..@ graphs      :List of 2
##   .. ..$ RNA_nn :Formal class 'Graph' [package "SeuratObject"] with 7 slots
##   .. .. .. ..@ assay.used: chr "RNA"
##   .. .. .. ..@ i         : int [1:454180] 0 995 1599 3985 8528 9237 9472 11547 11573 14618 ...
##   .. .. .. ..@ p         : int [1:22710] 0 12 42 67 71 73 91 117 141 152 ...
##   .. .. .. ..@ Dim       : int [1:2] 22709 22709
##   .. .. .. ..@ Dimnames  :List of 2
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. ..@ x         : num [1:454180] 1 1 1 1 1 1 1 1 1 1 ...
##   .. .. .. ..@ factors   : list()
##   .. ..$ RNA_snn:Formal class 'Graph' [package "SeuratObject"] with 7 slots
##   .. .. .. ..@ assay.used: chr "RNA"
##   .. .. .. ..@ i         : int [1:1509887] 0 995 1024 1156 1553 2310 3985 4234 5012 7282 ...
##   .. .. .. ..@ p         : int [1:22710] 0 40 106 192 235 291 348 432 494 595 ...
##   .. .. .. ..@ Dim       : int [1:2] 22709 22709
##   .. .. .. ..@ Dimnames  :List of 2
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. ..@ x         : num [1:1509887] 1 0.1429 0.1429 0.0811 0.1429 ...
##   .. .. .. ..@ factors   : list()
##   ..@ neighbors   : list()
##   ..@ reductions  :List of 3
##   .. ..$ pca    :Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
##   .. .. .. ..@ cell.embeddings           : num [1:22709, 1:50] 0.885 -1.466 -31.677 0.699 -1.37 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
##   .. .. .. ..@ feature.loadings          : num [1:2000, 1:50] -0.0236 -0.1313 -0.0414 0.0053 -0.1367 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:2000] "Chil3" "Cd74" "Lyz2" "Ccl5" ...
##   .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
##   .. .. .. ..@ feature.loadings.projected: num[0 , 0 ] 
##   .. .. .. ..@ assay.used                : chr "RNA"
##   .. .. .. ..@ global                    : logi FALSE
##   .. .. .. ..@ stdev                     : num [1:50] 6.34 5.34 4.72 4.46 4 ...
##   .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
##   .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ] 
##   .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ] 
##   .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ] 
##   .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ] 
##   .. .. .. ..@ misc                      :List of 1
##   .. .. .. .. ..$ total.variance: num 958
##   .. .. .. ..@ key                       : chr "PC_"
##   .. ..$ harmony:Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
##   .. .. .. ..@ cell.embeddings           : num [1:22709, 1:50] 2.142 0.244 -23.871 1.78 0.252 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. ..$ : chr [1:50] "harmony_1" "harmony_2" "harmony_3" "harmony_4" ...
##   .. .. .. ..@ feature.loadings          : num [1:2000, 1:50] -397 -966 -509 -1003 8498 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:2000] "Sox17" "St18" "Prex2" "Slco5a1" ...
##   .. .. .. .. .. ..$ : chr [1:50] "harmony_1" "harmony_2" "harmony_3" "harmony_4" ...
##   .. .. .. ..@ feature.loadings.projected: num [1:2000, 1:50] -397 -966 -509 -1003 8498 ...
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:2000] "Sox17" "St18" "Prex2" "Slco5a1" ...
##   .. .. .. .. .. ..$ : chr [1:50] "harmony_1" "harmony_2" "harmony_3" "harmony_4" ...
##   .. .. .. ..@ assay.used                : chr "RNA"
##   .. .. .. ..@ global                    : logi FALSE
##   .. .. .. ..@ stdev                     : num [1:50] 5.53 4.36 4.18 4.6 3.25 ...
##   .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
##   .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ] 
##   .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ] 
##   .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ] 
##   .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ] 
##   .. .. .. ..@ misc                      : list()
##   .. .. .. ..@ key                       : chr "harmony_"
##   .. ..$ umap   :Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
##   .. .. .. ..@ cell.embeddings           : num [1:22709, 1:2] 0.696 -4.174 12.588 -2.845 4.01 ...
##   .. .. .. .. ..- attr(*, "scaled:center")= num [1:2] -0.247 -0.207
##   .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. ..$ : chr [1:22709] "HDM6_6WK_AAACCTGCAATTCCTT-1" "HDM6_6WK_AAACCTGCACATTTCT-1" "HDM6_6WK_AAACCTGCATCCGGGT-1" "HDM6_6WK_AAACCTGGTACAGTGG-1" ...
##   .. .. .. .. .. ..$ : chr [1:2] "umap_1" "umap_2"
##   .. .. .. ..@ feature.loadings          : num[0 , 0 ] 
##   .. .. .. ..@ feature.loadings.projected: num[0 , 0 ] 
##   .. .. .. ..@ assay.used                : chr "RNA"
##   .. .. .. ..@ global                    : logi TRUE
##   .. .. .. ..@ stdev                     : num(0) 
##   .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
##   .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ] 
##   .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ] 
##   .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ] 
##   .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ] 
##   .. .. .. ..@ misc                      : list()
##   .. .. .. ..@ key                       : chr "umap_"
##   ..@ images      : list()
##   ..@ project.name: chr "HDMvsDA"
##   ..@ misc        : list()
##   ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
##   .. ..$ : int [1:3] 5 0 2
##   ..@ commands    :List of 8
##   .. ..$ NormalizeData.RNA             :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "NormalizeData.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:22"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "NormalizeData(object = HDMDA_Obj_filtered)"
##   .. .. .. ..@ params     :List of 5
##   .. .. .. .. ..$ assay               : chr "RNA"
##   .. .. .. .. ..$ normalization.method: chr "LogNormalize"
##   .. .. .. .. ..$ scale.factor        : num 10000
##   .. .. .. .. ..$ margin              : num 1
##   .. .. .. .. ..$ verbose             : logi TRUE
##   .. ..$ FindVariableFeatures.RNA      :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "FindVariableFeatures.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:23"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "FindVariableFeatures(.)"
##   .. .. .. ..@ params     :List of 12
##   .. .. .. .. ..$ assay              : chr "RNA"
##   .. .. .. .. ..$ selection.method   : chr "vst"
##   .. .. .. .. ..$ loess.span         : num 0.3
##   .. .. .. .. ..$ clip.max           : chr "auto"
##   .. .. .. .. ..$ mean.function      :function (mat, display_progress)  
##   .. .. .. .. ..$ dispersion.function:function (mat, display_progress)  
##   .. .. .. .. ..$ num.bin            : num 20
##   .. .. .. .. ..$ binning.method     : chr "equal_width"
##   .. .. .. .. ..$ nfeatures          : num 2000
##   .. .. .. .. ..$ mean.cutoff        : num [1:2] 0.1 8
##   .. .. .. .. ..$ dispersion.cutoff  : num [1:2] 1 Inf
##   .. .. .. .. ..$ verbose            : logi TRUE
##   .. ..$ ScaleData.RNA                 :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "ScaleData.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:24"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "ScaleData(.)"
##   .. .. .. ..@ params     :List of 9
##   .. .. .. .. ..$ assay             : chr "RNA"
##   .. .. .. .. ..$ model.use         : chr "linear"
##   .. .. .. .. ..$ use.umi           : logi FALSE
##   .. .. .. .. ..$ do.scale          : logi TRUE
##   .. .. .. .. ..$ do.center         : logi TRUE
##   .. .. .. .. ..$ scale.max         : num 10
##   .. .. .. .. ..$ block.size        : num 1000
##   .. .. .. .. ..$ min.cells.to.block: num 3000
##   .. .. .. .. ..$ verbose           : logi TRUE
##   .. ..$ RunPCA.RNA                    :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "RunPCA.RNA"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:35"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "RunPCA(.)"
##   .. .. .. ..@ params     :List of 10
##   .. .. .. .. ..$ assay          : chr "RNA"
##   .. .. .. .. ..$ npcs           : num 50
##   .. .. .. .. ..$ rev.pca        : logi FALSE
##   .. .. .. .. ..$ weight.by.var  : logi TRUE
##   .. .. .. .. ..$ verbose        : logi TRUE
##   .. .. .. .. ..$ ndims.print    : int [1:5] 1 2 3 4 5
##   .. .. .. .. ..$ nfeatures.print: num 30
##   .. .. .. .. ..$ reduction.name : chr "pca"
##   .. .. .. .. ..$ reduction.key  : chr "PC_"
##   .. .. .. .. ..$ seed.use       : num 42
##   .. ..$ Seurat..ProjectDim.RNA.harmony:Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "Seurat::ProjectDim.RNA.harmony"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:45"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr [1:2] "Seurat::ProjectDim(object, reduction = reduction.save, overwrite = TRUE, " "    verbose = FALSE)"
##   .. .. .. ..@ params     :List of 7
##   .. .. .. .. ..$ reduction      : chr "harmony"
##   .. .. .. .. ..$ assay          : chr "RNA"
##   .. .. .. .. ..$ dims.print     : int [1:5] 1 2 3 4 5
##   .. .. .. .. ..$ nfeatures.print: num 20
##   .. .. .. .. ..$ overwrite      : logi TRUE
##   .. .. .. .. ..$ do.center      : logi FALSE
##   .. .. .. .. ..$ verbose        : logi FALSE
##   .. ..$ FindNeighbors.RNA.harmony     :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "FindNeighbors.RNA.harmony"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:48"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "FindNeighbors(., reduction = \"harmony\", dims = 1:17)"
##   .. .. .. ..@ params     :List of 16
##   .. .. .. .. ..$ reduction      : chr "harmony"
##   .. .. .. .. ..$ dims           : int [1:17] 1 2 3 4 5 6 7 8 9 10 ...
##   .. .. .. .. ..$ assay          : chr "RNA"
##   .. .. .. .. ..$ k.param        : num 20
##   .. .. .. .. ..$ return.neighbor: logi FALSE
##   .. .. .. .. ..$ compute.SNN    : logi TRUE
##   .. .. .. .. ..$ prune.SNN      : num 0.0667
##   .. .. .. .. ..$ nn.method      : chr "annoy"
##   .. .. .. .. ..$ n.trees        : num 50
##   .. .. .. .. ..$ annoy.metric   : chr "euclidean"
##   .. .. .. .. ..$ nn.eps         : num 0
##   .. .. .. .. ..$ verbose        : logi TRUE
##   .. .. .. .. ..$ do.plot        : logi FALSE
##   .. .. .. .. ..$ graph.name     : chr [1:2] "RNA_nn" "RNA_snn"
##   .. .. .. .. ..$ l2.norm        : logi FALSE
##   .. .. .. .. ..$ cache.index    : logi FALSE
##   .. ..$ FindClusters                  :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "FindClusters"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:30:50"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "FindClusters(., resolution = 0.5)"
##   .. .. .. ..@ params     :List of 11
##   .. .. .. .. ..$ graph.name      : chr "RNA_snn"
##   .. .. .. .. ..$ cluster.name    : chr "RNA_snn_res.0.5"
##   .. .. .. .. ..$ modularity.fxn  : num 1
##   .. .. .. .. ..$ resolution      : num 0.5
##   .. .. .. .. ..$ method          : chr "matrix"
##   .. .. .. .. ..$ algorithm       : num 1
##   .. .. .. .. ..$ n.start         : num 10
##   .. .. .. .. ..$ n.iter          : num 10
##   .. .. .. .. ..$ random.seed     : num 0
##   .. .. .. .. ..$ group.singletons: logi TRUE
##   .. .. .. .. ..$ verbose         : logi TRUE
##   .. ..$ RunUMAP.RNA.harmony           :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
##   .. .. .. ..@ name       : chr "RunUMAP.RNA.harmony"
##   .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2024-08-02 16:31:03"
##   .. .. .. ..@ assay.used : chr "RNA"
##   .. .. .. ..@ call.string: chr "RunUMAP(., reduction = \"harmony\", dims = 1:17)"
##   .. .. .. ..@ params     :List of 25
##   .. .. .. .. ..$ dims                : int [1:17] 1 2 3 4 5 6 7 8 9 10 ...
##   .. .. .. .. ..$ reduction           : chr "harmony"
##   .. .. .. .. ..$ assay               : chr "RNA"
##   .. .. .. .. ..$ slot                : chr "data"
##   .. .. .. .. ..$ umap.method         : chr "uwot"
##   .. .. .. .. ..$ return.model        : logi FALSE
##   .. .. .. .. ..$ n.neighbors         : int 30
##   .. .. .. .. ..$ n.components        : int 2
##   .. .. .. .. ..$ metric              : chr "cosine"
##   .. .. .. .. ..$ learning.rate       : num 1
##   .. .. .. .. ..$ min.dist            : num 0.3
##   .. .. .. .. ..$ spread              : num 1
##   .. .. .. .. ..$ set.op.mix.ratio    : num 1
##   .. .. .. .. ..$ local.connectivity  : int 1
##   .. .. .. .. ..$ repulsion.strength  : num 1
##   .. .. .. .. ..$ negative.sample.rate: int 5
##   .. .. .. .. ..$ uwot.sgd            : logi FALSE
##   .. .. .. .. ..$ seed.use            : int 42
##   .. .. .. .. ..$ angular.rp.forest   : logi FALSE
##   .. .. .. .. ..$ densmap             : logi FALSE
##   .. .. .. .. ..$ dens.lambda         : num 2
##   .. .. .. .. ..$ dens.frac           : num 0.3
##   .. .. .. .. ..$ dens.var.shift      : num 0.1
##   .. .. .. .. ..$ verbose             : logi TRUE
##   .. .. .. .. ..$ reduction.name      : chr "umap"
##   ..@ tools       : list()
head(HDMDA_Obj_Harm@meta.data)
##                                orig.ident nCount_RNA nFeature_RNA   Id Time
## HDM6_6WK_AAACCTGCAATTCCTT-1 SeuratProject       1859          962 HDM6  6WK
## HDM6_6WK_AAACCTGCACATTTCT-1 SeuratProject       6871         2588 HDM6  6WK
## HDM6_6WK_AAACCTGCATCCGGGT-1 SeuratProject       3226         1227 HDM6  6WK
## HDM6_6WK_AAACCTGGTACAGTGG-1 SeuratProject       1637          842 HDM6  6WK
## HDM6_6WK_AAACCTGGTCAGATAA-1 SeuratProject       1919         1185 HDM6  6WK
## HDM6_6WK_AAACCTGGTCCAGTAT-1 SeuratProject       4150         1881 HDM6  6WK
##                                        Barcode PercentMT RNA_snn_res.0.5
## HDM6_6WK_AAACCTGCAATTCCTT-1 AAACCTGCAATTCCTT-1  2.850995               5
## HDM6_6WK_AAACCTGCACATTTCT-1 AAACCTGCACATTTCT-1  1.673701               0
## HDM6_6WK_AAACCTGCATCCGGGT-1 AAACCTGCATCCGGGT-1  4.463732              10
## HDM6_6WK_AAACCTGGTACAGTGG-1 AAACCTGGTACAGTGG-1  3.543067               4
## HDM6_6WK_AAACCTGGTCAGATAA-1 AAACCTGGTCAGATAA-1  2.397082               5
## HDM6_6WK_AAACCTGGTCCAGTAT-1 AAACCTGGTCCAGTAT-1  2.024096               6
##                             seurat_clusters
## HDM6_6WK_AAACCTGCAATTCCTT-1               5
## HDM6_6WK_AAACCTGCACATTTCT-1               0
## HDM6_6WK_AAACCTGCATCCGGGT-1              10
## HDM6_6WK_AAACCTGGTACAGTGG-1               4
## HDM6_6WK_AAACCTGGTCAGATAA-1               5
## HDM6_6WK_AAACCTGGTCCAGTAT-1               6
# Top 20 variable genes
top_Vgenes <- head(VariableFeatures(HDMDA_Obj_Harm), 20)
plot1 <- VariableFeaturePlot(HDMDA_Obj_Harm)+ theme(legend.position = "none")
plot2 <- LabelPoints(plot = plot1, points = top_Vgenes, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
plot1 + plot2
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.
## Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ElbowPlot(HDMDA_Obj_Harm)

PCHeatmap(HDMDA_Obj_Harm, dim = 1:6, cells = 500, balanced =TRUE, ncol = 3)

DimPlot(HDMDA_Obj_Harm, reduction = "umap", group.by = "RNA_snn_res.0.5", label = TRUE, label.box = TRUE)+ labs(title = "With Integration")

DimPlot(HDMDA_Obj_Harm, reduction = "umap", group.by = "RNA_snn_res.0.5", split.by = "Id", label = TRUE, label.box = TRUE)+ labs(title = "With Integration")

Annotate cell cluster

Cl_marker <- FindAllMarkers(HDMDA_Obj_Harm, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
## Calculating cluster 0
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 11
## Calculating cluster 12
## Calculating cluster 13
## Calculating cluster 14
## Calculating cluster 15
head(Cl_marker)
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster   gene
## Il1rl1     0   3.461304 0.814 0.127         0       0 Il1rl1
## Fgl2       0   2.290296 0.851 0.310         0       0   Fgl2
## Gata3      0   2.058888 0.884 0.383         0       0  Gata3
## Rbpj       0   1.473229 0.829 0.370         0       0   Rbpj
## Il13       0   5.526502 0.458 0.021         0       0   Il13
## Ltb4r1     0   2.074159 0.540 0.150         0       0 Ltb4r1
dim(Cl_marker)
## [1] 8359    7
#View top markers for the specific cluster
Top10 <- Cl_marker %>% group_by(cluster) %>% top_n(n=10, wt = avg_log2FC)
DoHeatmap(HDMDA_Obj_Harm, features = Top10$gene) + NoLegend()
## Warning in DoHeatmap(HDMDA_Obj_Harm, features = Top10$gene): The following
## features were omitted as they were not found in the scale.data slot for the RNA
## assay: Eaf2, Derl3, Fkbp11, Nedd4, Ttk, Espl1, Ankrd12, AC149090.1, Macf1,
## Neat1, AU020206, mt-Atp8, mt-Cytb, mt-Co3, Malat1, Klhdc1, B630019A10Rik,
## Asap1, Gm47242, Lpar6, 1700020I14Rik, Trp53i13, Hnrnpll, 2010111I01Rik, Nebl